Fastq-mcf -s parameter change

---

Project: Differences in pre-processing results if -s (scale) parameter of fastq-mcf is altered
Date: 08 November 2019
Report version: 0.2
Author(s): Kohl Kinning
kohl.kinning@cuanschutz.edu
Linda Crnic Institute for Down Syndrome, University of Colorado Anschutz Medical Campus

The -s (scale) parameter was set to 0 in the pipeline. This parameter is interpreted as log2(s), so the setting of 0 results in adapter trimming after a single base match. We noticed over-aggressive trimming of the ends of reads due to such a low threshold for matching. In this report, we initially look at the differences of several -s settings and then continue to look at the differences between a setting of 0 and 1.6.

Libraries

library(ggplot2)
library(ggpubr)
library(dplyr)
library(plotly)

Difference in total counts

A match of 1 and 2 bases looks similar, while higher settings result in a sharp drop off.

trim_counts <- read.table("~/Projects/CD4_CD8_NK_trimming_investigations/trimming_results.csv", header=TRUE, sep=",")

ggplot(trim_counts, aes(x=as.character(scale), y=clipped_end_reads)) +
  geom_bar(stat="identity") +
  xlab("-s parameter") + ylab("Clipped reads") +
  ggtitle("Effect of -s parameter on read clipping") +
  theme_pubr()

Read the counts data in

Look at the HTSeq resulting gene count data for the setting of 0 and 1.6.

zero_counts <- read.table("~/Projects/CD4_CD8_NK_trimming_investigations/counts/0/CD4_100B2_HTSeq_counts.txt", row.names=1)
one_six_counts <- read.table("~/Projects/CD4_CD8_NK_trimming_investigations/counts/1.6/CD4_100B2_HTSeq_counts.txt", row.names=1)

bound_counts <- zero_counts
names(bound_counts) <- "lower_thresh"
bound_counts$higher_thresh <- one_six_counts$V2
bound_counts$diff <- bound_counts$lower_thresh - bound_counts$higher_thresh
bound_counts$ratio_change <- bound_counts$diff/bound_counts$lower_thresh

Scatter plot

Observe no outliers when comparing the counts between the settings. Any changes are minor and in a consistent direction.

ggplotly(
ggplot(bound_counts, aes(x=log2(lower_thresh), y=log2(higher_thresh), label=rownames(bound_counts))) +
  geom_point(shape=1) +
  xlab("log2(-s 0 resulting counts)") + ylab("log2(-s 1.6 resulting counts)") +
  ggtitle("Difference in resulting counts\ndue to -s parameter") +
  theme_pubr()
)

Total difference in counts

Look at the total difference in HTSeq counts due to the -s parameter change

sum(bound_counts$lower_thresh)-sum(bound_counts$higher_thresh)

[1] 144020

Zero count genes

Number of genes lost completely

colSums(bound_counts[1:2]==0)

 lower_thresh higher_thresh 
         9305          9313 

Average change

Look at the average per-gene counts change, using the absolute differences

mean(abs(bound_counts$diff))

[1] 6.622256

Gene gain

Which genes do we gain with the higher thresh?

gained_genes <-
  setdiff(
    rownames(bound_counts[which(bound_counts$lower_thresh==0),]),
    rownames(bound_counts[which(bound_counts$higher_thresh==0),])
  )

bound_counts[gained_genes,1:3]

Gene loss

Which genes do we lose with the higher thresh?

lost_genes <-
  setdiff(
    rownames(bound_counts[which(bound_counts$higher_thresh==0),]),
    rownames(bound_counts[which(bound_counts$lower_thresh==0),])
  )

bound_counts[lost_genes,1:3]

Highest differences

Look at the genes which show the largest differences in counts